function [residuals,J] = nlsTheta1(theta1Values,setup)
[Q,output,errorMes] = objectFunc(theta1Values,setup);
if errorMes == 1
    select    = ~isnan(setup.yields);
    totalObs  = sum(sum(select));
    residuals = ones(totalObs,1)*1D60;
    numTheta1 = size(theta1Values,1);
    J = ones(totalObs,numTheta1)*1D-6;
   return
end
selectMat       = ~isnan(output.residuals);
totalObs        = sum(sum(selectMat));
% Scaling the results so residuals'*residuals correspnds to sqrt(sum(Qterms,1)/sum(nyIndex))*10
scalar          = 1/(sum(output.Qterms,1)*sum(output.nyIndex))^0.25*10;
residuals       = reshape(output.residuals(selectMat),totalObs,1)*scalar;

if nargout > 1
   % We compute the gradient with multiprocessing - using only 2 CPUs
   numTheta1 = size(theta1Values,1);
   J = ones(size(residuals,1),numTheta1)*1D-6;
   epsValue = max(abs(theta1Values)*1D-5,1D-10); 
   printToScreen(theta1Values)
   if setup.MultProcessOn == 1
       % We use multiprocessing
       parfor (i=1:numTheta1,2)
           theta1Values_p_eps       = theta1Values;
           theta1Values_p_eps(i,1)  = theta1Values(i,1) + epsValue(i,1);
           [~,output,errorMes]      = objectFunc(theta1Values_p_eps,setup);
           if errorMes == 0
               residuals_p           = reshape(output.residuals(selectMat),totalObs,1)*scalar;
               J(:,i)                = (residuals_p   - residuals)/epsValue(i,1);
           end
       end
   else
       % We use no multiprocessing
       for i=1:numTheta1
           theta1Values_p_eps       = theta1Values;
           theta1Values_p_eps(i,1)  = theta1Values(i,1) + epsValue(i,1);
           [~,output,errorMes]      = objectFunc(theta1Values_p_eps,setup);
           if errorMes == 0
               residuals_p           = reshape(output.residuals(selectMat),totalObs,1)*scalar;
               J(:,i)                = (residuals_p   - residuals)/epsValue(i,1);
           end
       end
   end
   % For elements of paramsValues close to the boundary, we downscale these columns
   select = (abs(setup.lowerBoundsValues-theta1Values) < 1D-5);
   J(:,select') = J(:,select')*1D-4;
   select = (abs(setup.upperBoundsValues-theta1Values) < 1D-5);
   J(:,select') = J(:,select')*1D-4;
end

end


